#pragma once
/******
 * Purpose: This program is design to do mapping SOLiD or Illumina Reads,
 * The input is a set of short reads file and a Genome_Indes_TableQ
 * The mapping result is output to a file
 *  Author: Yangho Chen
 */
#include "LongReadsSet.h"
#include "PairedReadsSet.h"
#include "ReadsMappingStats.h"
#include "MappingResult.h"
#include "Genome_Index_TableQ.h"
#include "ParameterList.h"
#include "AlignmentsQ.h"
#include "ReadInBitsSet.h"
#include "ReadInBits.h"
#include "ShortReadUtil.h"
#include "ColorSpaceRead.h"
#include "GenomeNTdata.h"
#include "chromosomeNTdata.h"
#include "stdafx.h"
#include <time.h>
#include <vector>
#ifdef WIN32
#include <direct.h>
#endif

// Use OpenMP is gcc version is later than 4.2
#ifdef __GNUC__
#ifdef __GNUC_PATCHLEVEL__
#define __GNUC_VERSION__ (__GNUC__ * 10000 \
                            + __GNUC_MINOR__ * 100 \
                            + __GNUC_PATCHLEVEL__)
#else
#define __GNUC_VERSION__ (__GNUC__ * 10000 \
                            + __GNUC_MINOR__ * 100)
# endif

#if __GNUC_VERSION__ >= 40200
#include <omp.h>
#endif
#else
#ifdef _MSC_VER
#if _MSC_VER >= 2000
#include <omp.h>
#endif
#endif
#endif

// Macro for Parallelization with OpenMP
#ifndef _OPENMP
#define __OPENMP_FOR_PARALLEL__(SHARP_SIGN, openmp_flag) \
{LOG_INFO("Info %d: Use single CPU because OpenMP is not available.\n", CONFIG_LOG);}
#else
#define __OPENMP_FOR_PARALLEL__(SHARP_SIGN, openmp_flag) {\
    int numberOfCPUs = omp_get_num_procs();\
    std::cout << numberOfCPUs << " CPUs." << BLANK_LINE << "\n";\
    SHARP_SIGN##openmp_flag omp parallel for\
}
#endif

// For both counter for chromosome and tolerated substitution error
const int ALIGNMENT_RESULT_FILE_BUFFER_SIZE = 1000000;

int parallelMappingLongReads(vector<string>& readSetsList,\
                             GIndexTQ& indexTable, MappingOpts P);

int parallelMapping(vector<string>& readSetsList,\
                    GIndexTQ& indexTable, MappingOpts P);

/*
 * This class maps Illumina or SOliD short reads to reference genome.
 * It switches different ways to do read mapping, according to different options.
 */
class CReadsMapping: public CReadsMappingStats
{
public:
    CReadsMapping(void);
    CReadsMapping(MappingOpts P);
    ~CReadsMapping(void);
    MappingOpts opt;
    int mapReadsSets(const char* ReadsSetsList, GIndexTQ& table);
    template <class ReadSetType>
    int mapReadsInAFile (ReadSetType& readSet, const GIndexTQ& table);
    int mapReads(CReadInBitsSet& readsSet, const GIndexTQ& table);
    int mapReads(CLongReadsSet& readsSet, const GIndexTQ& table);
    int queryARead(CReadInBits& r, const GIndexTQ& table, CAlignmentsQ& aQue) const;
    int queryALongRead(CReadInBits& r1stHalf, CReadInBits& r2ndHalf, const GIndexTQ& table, CAlignmentsQ& aQue) const;
protected:
    inline void getReadInfo(CReadInBitsSet& readSet, int readId, CReadInBits& r, CMappingResult& m);
    inline void getLongReadInfo(CLongReadsSet& readSet, int readId, \
                                CReadInBits& r1stHalf, CReadInBits& r2ndHalf, CMappingResult& m);
    int printMapInfo(CReadInBitsSet& readsSet, int seedOpt);
    int setUpIO4Aligment(const char* Solexafile, const GIndexTQ& table);
    int tearDownIO4Aligment(void);
    FileOutputBuffer*  AlignResult;
    FileOutputBuffer*  AmbiguousReads;
    FileOutputBuffer*  BadReads;
    FileOutputBuffer*  MissReads;
    time_t start, end;

    const static unsigned int CHECK_POINTS = 1000000;
    int dealMappedRead(const GIndexTQ& table, CAlignmentsQ& aQue);
    int dealMappedLongRead(const GIndexTQ& table, CAlignmentsQ& aQue, CMappingResult& m);
    int dealMissedRead(CMappingResult& m);
    int dealMissedRead(bool bMapReadInColors, const char* readName, CReadInBits r, const char* qs = NULL);
    int dealAmbiguousRead(CMappingResult& m);
    int dealAmbiguousRead(bool bMapReadInColors, const char* readName, CReadInBits r, const char* qs = NULL);
    int printSingleEndReads(CMappingResult& m);
    inline void printCheckPointInfo(int i);
    virtual int printLogFile(const char* inputFile);

    char cOutputFormat;
    inline bool isSamFormat() {
        return(this->cOutputFormat == 's');
    };
private:
    void initialization(void);
    inline void printRead(FileOutputBuffer* FileBuf, CMappingResult& m);
    inline void printRead(FileOutputBuffer* FileBuf, bool bMapReadInColors, const char* readName, CReadInBits r, const char* qs);
    inline void getLongBaseReadInfo(CLongReadsSet& readSet, int readId, CReadInBits& r1stHalf, CReadInBits& r2ndHalf, CMappingResult& m);
    inline void getLongColorReadInfo(CLongReadsSet& readSet, int readId, CReadInBits& r1stHalf, CReadInBits& r2ndHalf, CMappingResult& m);
    inline void getColorReadInfo(CReadInBitsSet& readSet, int readId, CReadInBits& r, CMappingResult& m);
    inline void getBaseReadInfo(CReadInBitsSet& readSet, int readId, CReadInBits& r, CMappingResult& m);
    string getMappingFileN(const char* caReadsSetName, const GIndexTQ& table);
    // unsigned int checkPairedReadSetSize(CReadInBitsSet& firstHalfSet, CReadInBitsSet& SecondHalfSet);
};

/*
 * Get read sequence and qualifty scores in m
 */
inline void CReadsMapping::getLongReadInfo\
(CLongReadsSet& readSet, int readId, CReadInBits& r1stHalf, CReadInBits& r2ndHalf, CMappingResult& m)
{
    if(this->opt.bMappedSOLiDRead) {
        getLongColorReadInfo(readSet, readId, r1stHalf, r2ndHalf, m);
    } else {
        getLongBaseReadInfo(readSet, readId, r1stHalf, r2ndHalf, m);
    }
};

inline void CReadsMapping::getReadInfo(CReadInBitsSet& readSet, int readId, CReadInBits& r, CMappingResult& m)
{
    if(this->opt.bMappedSOLiDRead) {
        getColorReadInfo(readSet, readId, r, m);
    } else {
        getBaseReadInfo(readSet, readId, r, m);
    }
}


inline void CReadsMapping::getLongBaseReadInfo\
(CLongReadsSet& readSet, int readId, CReadInBits& r1stHalf, CReadInBits& r2ndHalf, CMappingResult& m)
{
    CReadInBitsSet& readSet1stHalf = *(readSet.F_Reads);
    CReadInBitsSet& readSet2ndHalf = *(readSet.R_Reads);
    m.uiReadLength = readSet.getReadLength(); // Check the read length
    // Get read tag
    if (readSet.size() > readId) {
        strcpy(m.QNAME, readSet1stHalf.pReadsID->at(readId).id);
    } else {
        sprintf(m.QNAME, "Read_%d", this->iReadCounter + readId);
    }
    // Get read sequence
    // TODO: modify this->opt because the anchorLength could be anchorLength2
    unsigned int halfReadLength = this->opt.anchorLength;
    unsigned int secondHalfStart = readSet.getReadLength() - halfReadLength;
    r1stHalf.decode(m.caRead);
    r2ndHalf.decode(&(m.caRead[secondHalfStart]));
    // Get the quality score
    if (readSet1stHalf.pQualScores != NULL) {
        const char* qS1stHalf = readSet1stHalf.pQualScores->qScores((unsigned int)readId);
        trQScores(halfReadLength, 0, qS1stHalf, m.rawScores);
        const char* qS2ndHalf = readSet2ndHalf.pQualScores->qScores((unsigned int)readId);
        trQScores(halfReadLength, 0, qS2ndHalf, &m.rawScores[secondHalfStart]);
        trQScores(m.uiReadLength, SolexaScoreEncodingShift, m.rawScores, m.QScores);
    }
    if (this->cOutputFormat == 's') {
        m.getReverseReadandQual();
    }
}

inline void CReadsMapping::getLongColorReadInfo\
(CLongReadsSet& readSet, int readId, CReadInBits& r1stHalf, CReadInBits& r2ndHalf, CMappingResult& m)
{
    CReadInBitsSet& readSet1stHalf = *(readSet.F_Reads);
    CReadInBitsSet& readSet2ndHalf = *(readSet.R_Reads);
    m.uiReadLength = readSet.getReadLength(); // Check the read length
    // Get read tag
    if (readSet.size() > readId) {
        strcpy(m.QNAME, readSet1stHalf.pReadsID->at(readId).id);
    } else {
        sprintf(m.QNAME, "Read_%d", this->iReadCounter + readId);
    }
    // Get read sequence
    // TODO: the anchorLength may be anchorLength2 if the 2nd end is also long
    unsigned int halfReadLength = this->opt.anchorLength;
    unsigned int secondHalfStart = readSet.getReadLength() - halfReadLength;
    // TODO: decodeColors is to get colors but the base should be the one precent in buffer?
    CLongReadInBits r(r1stHalf, r2ndHalf);
    decodeColors(m.caRead, r, readSet.getReadLength());
    // Get the quality score for SOLiD
    /* TODO fix the SOLiD score (base score and color space score*/
    if (readSet1stHalf.pQualScores != NULL) {
        const char* qS1stHalf = readSet1stHalf.pQualScores->qScores((unsigned int)readId);
        trQScores(halfReadLength, 0, qS1stHalf, m.rawScores);
        const char* qS2ndHalf = readSet2ndHalf.pQualScores->qScores((unsigned int)readId);
        trQScores(halfReadLength, 0, qS2ndHalf, &m.rawScores[secondHalfStart]);
        trQScores(m.uiReadLength, Phred_SCALE_QUAL_SHIFT, m.rawScores, m.QScores);
    }
    // For sam Format, one need reversed quality and Seq
    if (this->cOutputFormat == 's') {
        m.getReverseReadandQual();
    }
}

inline void CReadsMapping::getBaseReadInfo(CReadInBitsSet& readSet, int readId, CReadInBits& r, CMappingResult& m)
{
    m.uiReadLength = readSet.getReadLength(); // Check the read length
    // Get read tag
    if ((int)readSet.pReadsID->size() > readId) {
        strcpy(m.QNAME, readSet.pReadsID->at(readId).id);
    } else {
        sprintf(m.QNAME, "Read_%d", this->iReadCounter + readId);
    }
    // TODO: Modified r.decode depends on the length of CReadInBits.
    r.decode(m.caRead);
    // Get the quality score for SOLEXA
    if (readSet.pQualScores != NULL) {
        const char* qS = readSet.pQualScores->qScores((unsigned int)readId);
        trQScores(readSet.getReadLength(), 0, qS, m.rawScores);
        trQScores(readSet.getReadLength(), SolexaScoreEncodingShift, m.rawScores, m.QScores);
    }
    if (this->cOutputFormat == 's') {
        m.getReverseReadandQual();
    }
}

inline void CReadsMapping::getColorReadInfo(CReadInBitsSet& readSet, int readId, CReadInBits& r, CMappingResult& m)
{
    m.uiReadLength = readSet.getReadLength(); // Check the read length
    // Get read tag
    if ((int)readSet.pReadsID->size() > readId) {
        strcpy(m.QNAME, readSet.pReadsID->at(readId).id);
    } else {
        sprintf(m.QNAME, "Read_%d", this->iReadCounter + readId);
    }
    decodeColors(m.caRead, r);
    // Get the quality score for SOLiD
    /* TODO fix the SOLiD score (base score and color space score*/
    if (readSet.pQualScores != NULL) {
        const char* qS = readSet.pQualScores->qScores((unsigned int)readId);
        trQScores(readSet.getReadLength(), 0, qS, m.rawScores);
        trQScores(readSet.getReadLength(), Phred_SCALE_QUAL_SHIFT, m.rawScores, m.QScores);
    }
    if (this->cOutputFormat == 's') {
        // For sam Format, one need reversed quality and Seq
        m.getReverseReadandQual();
    }
}

inline void CReadsMapping::printCheckPointInfo(int readNo)
{
    if (readNo % this->CHECK_POINTS == 0) {
        printf("Mapping no %u reads.\r", this->iReadCounter + readNo);
        fflush(stdout);
    }
}

// The following function fill different part of CMappingResult
inline
bool getSingleMappingIndex(CGenomeNTdata& pgenomeNT, CAlignmentsQ& aQue, int mappingIndex, CMappingResult &m)
{
    m.uiDiff = aQue.asdiff[mappingIndex];
    m.MultipleMappedNo = aQue.load;
    m.strand = (mappingIndex >= (int)aQue.ForwardAlignmentLoad) ? '-' : '+' ;

    m.uiGlobalMappedPos = aQue.aiHitIndex[mappingIndex];
    m.uiRefId = pgenomeNT.genomeIndex2chrID(m.uiGlobalMappedPos);
    m.uiPOS = pgenomeNT.genomeLocusID2chrIndex(m.uiGlobalMappedPos);
    ChrIndex2GeneName& geneVec = pgenomeNT.paChromosomes[m.uiRefId]->geneVec;
    if (geneVec.table.size() > 0) {
        CGene g = geneVec.query(m.uiPOS);
        strcpy(m.RNAME, g.name.c_str());
        if(!g.isValid) { // this is for handeling mapping exception for mapping 2 NULL_REGION
            char* chrName = pgenomeNT.paChromosomes[m.uiRefId]->caInputFileName;
            char* strBuf = &m.RNAME[strlen(g.name.c_str())];
            myStrCpy(strBuf, chrName, FILENAME_MAX/2);
            return(false);
        }
        // g.startIndex is the translated index of m.uiPOS. Not the start index of gene
        m.uiPOS = g.startIndex;
    } else {
        sprintf(m.RNAME, "%d", m.uiRefId);
    }
    return(true);
}
void getQscores4Solexa(CAlignmentsQ& aQue, CMappingResult& m, bool samFormat);
void getReadQscores4Solexa(CAlignmentsQ& aQue, CMappingResult& m, bool samFormat);
void getSingleMappingSeq4Solexa(const GIndexTQ& table, CMappingResult& m, bool samFormat);
void getLongMappingInfo(const GIndexTQ& table, CAlignmentsQ& aQue, unsigned int mappingId, CMappingResult& m, \
                        bool samFormat, bool bNH);
void getSingleMappingInfo(const GIndexTQ& table, CAlignmentsQ& aQue, unsigned int mappingId, CMappingResult& m,\
                          bool samFormat, bool bNH);

inline void printSamHeader(FileOutputBuffer* AlignResult, vector<CGene>&refs, const char* RG, const char* CL)
{
    sprintf(AlignResult->caBufp, "@HD\tVN:0.1.5\tSO:queryname\n");
    AlignResult->UpdateSize();
    for (vector<CGene>::iterator it = refs.begin(); it != refs.end(); it++ ) {
        sprintf(AlignResult->caBufp, "@SQ\tSN:%s\tLN:%u\n", it->name.c_str(), it->startIndex);
        AlignResult->UpdateSize(); // startIndex is actually the length of reference.
    }
    sprintf(AlignResult->caBufp, "%s\n@PG\tID:PerM\tVN:0.4.1\tCL:\"%s\"\n", RG, CL);
    AlignResult->UpdateSize();
}

inline const char* getLongRefSeq(const GIndexTQ& table, CMappingResult& m, bool bNoRef)
{
    if(bNoRef) {
        m.caRef[0] = '\0';
    } else {
        // TODO: Check and rename table.uiRead_Length
        unsigned int secondHalfStart = m.uiReadLength - table.uiRead_Length;
        CReadInBits ref1stHalf = table.pgenomeNTInBits->getSubstringInBits\
                                 (m.uiGlobalMappedPos, table.uiRead_Length);
        CReadInBits ref2ndHalf = table.pgenomeNTInBits->getSubstringInBits\
                                 (m.uiGlobalMappedPos + secondHalfStart, table.uiRead_Length);
        ref1stHalf.decode(m.caRef);
        ref2ndHalf.decode(&(m.caRef[secondHalfStart]));
        if(m.strand == '-') {
            reverseComplementKmer(m.caRef); // reverse complement reference
        }
    }
    return(m.caRef);
}
template <class ReadSetType>
bool wrongIndex(ReadSetType& readsSet, const GIndexTQ& table)
{
    // TODO the following check doesn't work for CLongReadsSet
    /*
    if (readSet.uiRead_Length != table.uiRead_Length) {
        LOG_INFO("Info %d: The index is not for read length %d, is for %d.\n",\
                 readSet.uiRead_Length, table.uiRead_Length, ERROR_LOG);
        return(true);
    }*/
    // TODO: check why .isSOLiDReadFormat() doesn't work and if cFileType is set correctly.
    bool bColorRead = readsSet.isSOLiDReadFormat();
    if ( bColorRead != table.bMapReadInColors) {
        if (bColorRead) {
            LOG_INFO("Info %d: The index is not for color read.\n", ERROR_LOG);
        } else {
            LOG_INFO("Info %d: The index is for color read.\n", ERROR_LOG);
        }
        return(true);
    }
    return(false);
}

template <class ReadSetType1, class ReadSetType2>
unsigned int checkPairedReadSetSize\
(ReadSetType1& firstHalfSet, ReadSetType2& SecondHalfSet)
{
    unsigned int size1 = (unsigned int)firstHalfSet.size();
    unsigned int size2 = (unsigned int)SecondHalfSet.size();
    if (size1 == size2) {
        return(size1);
    } else {
        LOG_INFO("Info %d: Not every read has the second half\n", WARNING_LOG);
        return(min(size1, size2));
    }
}

